How glassy are orientational dynamics of rodlike molecules near the isotropic-nematic 

transition? 
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In an attempt to quantitatively characterize the recently observed slow dynamics in the isotropic 
and nematic phase of liquid crystals, we investigate the single-particle orientational dynamics of 
rodlike molecules across the isotropic-nematic transition in computer simulations of a family of 
model systems of thermotropic liquid crystals. Several remarkable features of glassy dynamics 
are on display including non-exponential relaxation, dynamical heterogeneity, and non-Arrhenius 
temperature dependence of the orientational relaxation time. In order to obtain a quantitative 
measure of glassy dynamics in line with the estbalished methods in supercooled liquids, we construct 
a relaxation time versus scaled inverse temperature plot, and demonstrate that one can indeed define 
a 'fragility index' for thermotropic liquid crystals, that depends on density and aspect ratio. The 
values of the fragility parameter are surprisingly in the range one observed for glass forming liquids. 
A plausible correlation between the energy landscape features and the observed fragility is discussed. 

PACS numbers: 61.20.Lc,61.30.-v,64.70.Pf 



INTRODUCTION 



Thermotropic liquid crystals exhibit exotic phase behavior upon temperature variation. The nematic phase is rich 
\ with a long-ranged orientational order but lacks translational order. The isotropic-nematic (TN) phase transition, 
£5 ■ which is believed to be weakly first order in nature with certain characteristics of the continuous transition, has been 
a subject of immense attention in condensed matter physics and material sciences [H, @|- I n contrast, the dynamics 
of thermotropic liquid crystals have been much less studied, the focus being mostly on the long-time behavior of 
orientational relaxation near the TN transition A series of optical Kerr effect (OKE) measurements have, however, 
recently studied collective orientational relaxation in the isotropic phase near the TN transition over a wide range of 
time scales Q. The dynamics have been found to be surprisingly rich, the most intriguing feature being the power law 
decay of the OKE signal at short-to-intermediate times [3| • The relaxation scenario appears to be strikingly similar to 



that of supercooled molecular liquids Q , even though the latter do not undergo any thermodynamic phase transition. 
CO ■ Although the analogous dynanucs haVe been investigated in subsequent studies &|, a quantitative estimation of 
<0 glassy dynamics of rodlike molecules near the TN transition still eludes us. 



The prime objective of this paper is to provide a quantitative measure of glassy dynamics near the I-N transition. 
To this end, we have undertaken molecular dynamics simulations of a family of model systems consisting of rodlikc 
molecules across the TN transition in search of glassy behavior. Given the involvement of the phase transition to an 
orientationally ordered mesophase upon lowering the temperature, we choose to probe the single-particle orientational 
I ■ dynamics. We have defined a 'fragility index' and explored plausible correlation of the features of the underlying energy 
' O ' landscape with the observed fragility in analogy with supercooled liquids. 
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MODEL AND DETAILS OF THE SIMULATION 



The systems we have studied consist of ellipsoids of revolution. The Gay-Berne (GB) pair potential [7|, that is well 
established to serve as a model potential for systems of thermotropic liquid crystals, has been employed. The GB pair 
potential, which uses a single-site representation for each ellipsoid of revolution, is an elegant generalization of the 
extensively used isotropic Lennard- Jones potential to incorporate anisotropy in both the attractive and the repulsive 
parts of the interaction 0, Q.In the GB pair potential, ith ellipsoid of revolution is represented by the position of 
its center of mass and a unit vector e, along the long axis of the ellipsiod. The interaction potential between two 
ellipsoids of revolution i and j is given by 
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Here a ss defines the thickness or equivalently, the separation between the two ellipsoids of revolution in a side-by-side 



configuration, r.y is the distance between the centers of mass of the ellipsoids of revolution i and j , and 



i/n 



is a unit vector along the intermolecular separation vector r^. The molecular shape parameter a and the energy 
parameter e both depend on the unit vectors and as well as on fy as given by the following set of equations: 



a(rij,ei,ej) = <tq 
with x = (k 2 - 1)/(k 2 + 1) and 
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where the exponents p and ^ are adjustable parameter, and 
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withx' = W 1 ' i1 -1)/{k! + Here k = <j ee /a ss is the aspect ratio of the ellipsoid of revolution with a ee denoting 
the separation between two ellipsoids of revolution in a end-to-end configuration, and a ss = a$, and k' — e ss /e ee , 
where e ss is the depth of the minimum of the potential for a pair of ellipsoids of revolution aligned in a side-by-side 
configuration, and e ee is the corresponding depth for the end-to-end alignment, eo is the depth of the minimum of 
the pair potential between two ellipsoids of revolution alligncd in cross configuration. The GB pair potential defines 
a family of models, each member of which is characterized by the values chosen for the set of four parameters k, k', /i, 
and v, and is represented by GB(k,k',/j, ^) [8j. Three systems, namely GB(3, 5, 2, 1), GB(3.4, 5, 2, 1), and GB(3.8, 
5, 2, 1), that differ in the aspect ratio have been investigated. Molecular dynamics simulations have been performed 
with each of these systems, consisting of 500 ellipsoids of revolution, in a cubic box with periodic boundary conditions. 
Each of these systems has been studied along three isochors (p — 0.31, 0.32, and 0.33 for k — 3.0; p — 0.25, 0.26, and 
0.27 for k = 3.4; p = 0.215, 0.225, and 0.235 for k — 3.8) at several temperatures, starting from the high-temperature 
isotropic phase down to the nematic phase across the I-N phase boundary. All quantities are given in reduced units 
defined in terms of the Gay-Berne potential parameters eo and ctq: length in units of cto, temperature in units of 

and time in units of (^p) 1 ^ 2 , m being the mass of the ellipsoids of revolution. The mass as well as the moment of 
inertia of each of the ellipsoids of revolution have been set equal to unity. The intermolecular potential is truncated 
at a distance r cu tand shifted such that U(rij = r cu t) = 0, being the separation between two ellipsoids of revolution 
i and j . The equations of motion have been integrated using the velocity- verlet algorithm with integration time step 
dt = 0.0015 (9(. Equilibration has been done by periodic rescaling of linear and angular velocities of particles. This 
has been done for a time period of t q following which the system has been allowed to propagate with a constant 
energy for a time period of t e in order to ensure equilibration upon observation of no drift of temperature, pressure, 
and potential energy. The data collection has been executed in a microcanonical ensemble. At each state point, local 
potential energy minimization has been executed by the conjugate gradient technique for a subset of 200 statistically 
independent configurations. The landscape analysis has been done with a system size of 256 ellipsoids of revolution, 
which is big enough for having no qualitative change in the results due to the system size [l Oj ] . Minimization has 
been performed with three position coordinates and two Euler angles for each particle, the third Euler angle being 
redundant for ellipsoids of revolution. 



RESULTS AND DISCUSSION 



The single-particle second rank orientational time correlation function (OTCF) C|(t) is defined by 

<Ei^2(e i (t)-e i (0) > 



cm = 
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where P2 is the second rank Legendre polynomial, is the unit vector along the long axis of ith ellipsoid of revolution, 
and the angular brackets stand for ensemble averaging. 



FIG. 1: Time evolution of the single-particle second rank orientational time correlation function CJ(t) in a log-log plot for the 
system with aspect ratio k — 3. The time dependence is shown at several temperatures across the isotropic-nematic transition 
along the isochor at density p = 0.33. The solid lines denote the curves for the isotropic phase and the dashed lines for the 
nematic phase. The inset shows the average orientational order parameter < S > as a function of temperature. 



Figure Q] shows the time evolution of the single-particle second rank OTCF for one of the three systems considered 
here as the temperature is lowered along an isochor from the high-temperature isotropic phase down to the nematic 
phase across the FN phase boundary. In the inset, the average orientational order parameter < S > is shown as a 



function of temperature along the isochor The variation of < S > with temperature serves to locate the FN 



phase boundary. In the present study, the FN transition temperature Tj_jv is taken as the temperature at which 
< S > of the system is 0.35. For each aspect ratio, three isochors at different densities have been considered. The 
qualitative behavior has been found to be the same for all the three systems along all the isochors studied (data not 
shown). The emergence of the power law decay in the isotropic phase near the FN transition is evident in all the 
cases as a universal characteristic of FN transition [lj|. As the FN phase boundary is crossed upon cooling, the 
advent of two power law decay regimes separated by an intervening plateau at short-to-intermediate times imparts a 
step-like feature to the temporal behavior of the second rank OTCF. Such a feature bears remarkable similarity to 
what is observed for supercooled liquids as the glass transition is approached from the above 13j, LLJ] . While for the 
supercooled liquid the emargence of step-like feature in the OTCF is well understood as a consequence of relaxation, 
the origin of such a feature observed for liquid crystal defied of reliable explanation. 

We estimate the orientational correlation time r as the time taken for C|(t) to decay by 90%, i.e., Cf (i = r) = 0.1. 
Figure Ufa) shows t in the logarithmic scale as a function of the inverse temperature along the three isochors for each 
of the three systems considered. We have scaled the temperature by Ti-n in the spirit of Angell's plot, that displays 
the shear viscosity (or the structural relaxation time, the inverse diffusivity, etc.) of glass-forming liquids as a function 
of the inverse of the scaled temperature, the scaling being done in the latter case by the glass transition temperature 
T g 15, [l6|. For all the three systems, two distinct features are common: (i) in the isotropic phase far away from 



the FN transition, the orientational correlation time r exhibits the Arrhenius temperature dependence, i.e., t(T) 
Toexp[E / (ksT)], where the activation energy E and the pre-factor tq are both independent of temperature; (ii) in the 
isotropic phase near the FN transition, the temperature dependence of r shows marked deviation from the Arrhenius 
behavior and can be well described by the Vogel-Fulcher-Tammann (VFT) equation t(T) = r exp[B/(T — T VFT )], 
where tq, B, and Tvft are constants, independent of temperature. Again these features bear remarkable similarity 
with those observed for fragile glass-forming liquid. A non- Arrhenius temperature behavior is taken to be the signature 
of fragile liquids. For fragile liquids, the temperature dependence of the shear viscosity follows the Arrhenius behavior 
far above T g and can be fitted to the VFT functional form in the deeply supercooled regime near T g HI, 



The striking resemblance in the dynamical behavior described above between the isotropic phase of thermotropic 
liquid crystals near the FN transition and supercooled liquids near the glass transition has prompted us to attempt a 
quantitative measure of glassy behavior near the FN transition. For supercolled liquids, one quantifies the dynamics 
by a parameter called fragility index which measures the rapidity at which the liquid's propeties (such as viscosity) 
change as the glassy state is approached. In the same spirit 17j that offers a quantitative estimation of the fragile 
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FIG. 2: (a) The orientational correlation time r in the logarithmic scale as function of the inverse of the scaled temperature, the 
scaling being done by the isotropic to nematic transition temperature Tj_jv. For the insets, the horizontal and the vertical axis 
labels read same as that of the main frame and are thus omitted for clarity. Along each isochor, the solid line is the Arrhenius 
fit to the subset of the high-temperature data and the dotted line corresponds to the fit to the data near the isotropic-nematic 
phase boundary with the VFT form, (b) The fragility index m shown as a function of density for different aspect ratios. The 
dashed lines are guide to the eye to illustrate that the change in the fragility index for a given density difference (Ap) increases 
with the decrease in the aspect ratio. 



behavior of supercooled liquids, we here define the fragility index m of a thermotropic liquid crystalline system as 

dlogr(T) 



rn 



d(T!_ N /T) 



(8) 



T=T, 



Figure EJb) shows the density dependence of the fragility index for the three systems with different aspect ratios. For 
a given aspect ratio, the fragility index increases with increasing density, the numerical values of the fragility index 
m being comparable to those of supercooled liquids. The change in the fragility index for a given density difference 
(Ap) increases with the decrease in the aspect ratio. 



Another hallmark of fragile glass-forming liquids is spatially heterogeneous dynamics 18J reflected in non-Gaussian 
dynamical behavior [l9j |. It is intuitive that the growth of the pseudo- nematic domains, characterized by local nematic 
order, in the isotropic phase near the I-N transition would result in heterogeneous dynamics in liquid crystals. We 
have, therefore, monitored the time evolution of the rotational non-Gaussian parameter (NGP) [20(, a^{t), which in 
the present case is defined as 

o < A^ 4 (t) > -i m 



2 < A0 2 (t) > 5 



where 



< A^"(t) >= i J2 < \Mt) M0)\ 2n > • (io) 



Here 4>i is the rotation vector like the position vector appears incase of translational NGP of ith ellipsoid of 



revolution, the change of which is defined by A0i(t) = <j>iif) — </>i(0) = J Q dt'u(t'), uji being the corresponding angular 



velocity [13|, [14|, and N is the number of ellipsoids of revolution in the system. NGP will have value equal to zero 
when system dynamics is spatially homogeneous and will have a non-zero value when the system dynamics is spatially 
heterogeneous. As a typical behavior, Fig. EJa) shows the time dependence of the rotational NGP for one of the 
systems at several temperatures across the I-N transition along an isochor. On approaching the I-N transition upon 
cooling, a bimodal feature starts appearing with the growth of a second peak, which eventually becomes the dominant 
one, at longer times. 

We further investgate the apearance of this bimodal feature in NGP plot. To this end we calculate mean square 
angular deviation (MSAD) of the system at different temperatures starting from high temperature isotropic phase to 
low temperature nematic phase. The appearance of the bimodal feature in the rotational NGP is accompanied by 
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FIG. 3: (a) Time evolution of the rotational non-Gaussian parameter a 2 (t) in a semi-log plot for the system with aspect ratio 
k = 3. The time dependence is shown at several temperatures across the isotropic-nematic (I-N) transition along an isochor 
at density p — 0.33. (b) On a different scale along the vertical axis (appearing on the right), time evolution of the mean 
square angular deviation < A<p 2 (t) > is shown in a log- log plot for three temperatures: the highest temperature studied in 
the isotropic phase and the other two temperatures that are nearest to the I-N transition in the from both side along with the 
time evolution of a^it), and (c) On a different scale along the vertical axis (appearing on the right), the time evolution of the 
single-particle second rank orientational time correlation function CJ(t) is shown in a log-log plot for the two temperatures 
that are closest to the I-N transition on either side along with the time evolution of a^it). 



a signature of a sub-diffusive regime in the temporal evolution of the mean square angular deviation (MS AD), the 
time scale of the short-time peak and that of the onset of the sub-diffusive regime being comparable, as shown in 
Fig. [3Jb). We further note that the dominant peak appears on a time scale which is comparable to that of onset of 
the diffusive motion in orientational degrees of freedom (ODOF) as evident in Fig. EJb) and similar feature has been 
observed recently for glassy systems [20( also. 

We also observe that the time scale at which long-time peak appears is also comparable to the onset of the plateau 
that is observed in the time evolution of C|(t), as shown in Fig. [3Jc). 



Mode coupling theory analysis 

This striking similarities of the dynamics between liquid crystals near I-N transition and supercooled liquid near 
glass transition are also supported by the use of Mode Coupling Theory (MCT) to explain the dynamics of both the 
systems. While MCT was used first for the supercooled liquid, recently it has been used for liquid crystals also. MCT 
theory devloped by Gottke et al. @ predicts that near I-N transition, the low frequency rotational memory kernel 
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should diverge in a power law fashion. 



M R {z) « 4 ("J 



Mean field treatment gives a = 0.5. Invoking the rank (I) dependence of the memory function, the single particle 
OTCF can be written as [HI, HI H3 



Cf{z) = 



l(l + l)k B T]' ] 

J[* + r,(«) 



(12) 



The above equation can be Laplace inverted to obtain a short-to-intermadiate power law decay in C|(i) which is a 
universal characteristic of the I-N transition for several model liquid crystals [H, [(1 EH ■ Recently, Li et al. [25] have 
showed that it is also possible to formulate a schematic model that combines short-to-intermediate time relaxation 
with long time relaxation. In their model, they have expressed the total memory functon (M(t)) as the sum of mode 
coupling memory function (MMCT(t)) and Landau-de Gennes memory function (MtdGf(i))- 

M(t) = M MCT (t) + M LdG {t) (13) 

where 

M MCT (t) = n 2 K(t) (14) 

Q is the characteristic frequecny and K(0) = 1. Time dependence of K(t) can be written expressed in terms of the 
memory function m(t) [25} and m(t) has the following form 

m(t) = K<t>(t)<i>i(t) (15) 

Where <j>(i) is the autocorrelation function of the anisotropy of the polarizability and <j>i(t) is the solution of a F12 
schematic model for what is referred to as the density correlator, k being the coupling constant between them. Now, 
MLdc(t) can be written as 

M LdG (t) = TS(t) (16) 
Here is the relaxation time (r _1 = Tl<ig) and it diverges as (T — T*) _1 as the critical temperature T* of the I-N 



transition is approached from the above. Following the calculation of Ref. 25j , one can get two important relaxation 
equattions. 



-.2 



&(t) = -n;^i(t) -mMt) -n? / dt' mi (t-t')Mt') (17) 

Jo 

With the initial conditions <\>\ (0) = 1 and <j>\ (0) = and 

<j>(t) = -{n 2 + nr)(f>(t) - (n + r)<j>(t) - n 2 [ dt'm(t - t')4>{t') - n 2 r [ dt'm{t - t')<j)(t') (is) 

Jo Jo 

with the initial conditions 0(0) = 1 and 0(0) = — T . /ii and /j, are the damping constants. Eq. 17 is identical to one 
gets from the MCT analysis of the supercooled liquids. Difference between this scematic model and one applied for 
the supercooled liquid is in Eq. 18 . If T is set equal to in Eq. 18 , the supercooled model is recovered. Eq. 18 is 
the orientational correlation function coupled to the density correlation function with specific new terms that account 
for the long-time portion of the relaxation profile that has been previously described by LdG theory. 



Energy landscape analysis 



Several studies have attempted to interpret the fragility of glass-forming liquids in terms of the features of the 
underlying energy landscapes 2(| 27 , 28 , H^, [3(J 31 1 . Energy landscape analysis gives the potential energy, which 



devoids of any kind of thermal motions, of inherent structures of the parent liquid and hence provide a better 
understanding of the structure and dynamics of the paprent liquid. A recent study on thermotropic liquid crystals 
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FIG. 4: The energy landscapes as explored by the system with the aspect ratio 3 on variation of temperature along three 
isochors. (a) The average inherent structure (IS) energy per particle as a function of the inverse temperature at different 
densities. The solid line, dashed line, and dot-dashed lines are the linear fits to the data at the densities p — 0.31,0.32,0.33, 
respectively, (b) The displaced average IS energy per particle versus a scaled inverse temperature along the same three isochors. 
If the Gaussian model for the number density of the IS energy with a given energy is validated, a collapse of the data for all 
densities is expected onto a straight line with negative unit slope, that is drawn. The data are shown over the temperature 
regimes in which the average IS energy is on a decline. 



has reported the temperature dependent exploration of the energy landscapes of a family of the Gay-Berne model 
systems across the mesophases [10|. The average inherent structure (IS) energy < e/s > has been found to fall as 
< S > grows across the I-N phase boundary and through the nematic phase in contrast to its insensitivity to the 
temperature in the high-temperature isotropic phase and the low-temperature smectic-B phase [10]. Such a fall in the 
average IS energy is consistent with a Gaussian form for the number density of inherent structures with energy ejs, 
that predicts a linear variation of < ejs > with the inverse temperature: < eis > (T) = e° T g — a 2 /2NksT, where 
e° IS and a are parameters independent of temperature and ks is the Boltzmann constant [29(. Note that this has 
been observed for a glassy system [29j . where the average IS energy also falls over a temperature range [27| • In Fig. 
H]Ja), we demonstrate with the original and the most studied parameterization for the GB pair potential GB(3, 5, 2, 1) 
that the prediction holds good over the temperature range where < eis > is on a decline along all the three isochors 
studied. It then follows that a plot of < eis > {T) — e° IS versus a 2 /2NksT would result in a collapse of the < eis > 
data for all densities onto a straight line with negative unit slope. This is indeed found to be true, as shown in Fig. 
HJb), implying the validity of the Gaussian model in this case as well. It may be noted that when the distribution of 
IS energy is Gaussian, the fragility of glass-forming liquids has been shown to depend on the total number of inherent 
structures, the width of the Gaussian, and the variation of the basin shape with the average IS energy [291 ]. 



CONCLUSION 



The origin of the glassy orientational dynamics, both single particle and collective, of nematogens near the isotropic- 
nematic transition has been addressed in several publications in recent years [3|, |4| , la, llQj • In these studies, the similarity 
between dynamics of supercooled glassy liquid and liquid crystals has been discussed in detail, but no quantitative 
measure of the similarity was provided. The fragility index introduced here serves to remove that lacuna. It is indeed 
surprising that even the values of the fragility parameter are in the range observed for glassy liquids as well. This is in 
agreement with the repeated observation by Fayer and coworkers that the values of the power law exponents observed 
in the two systems are quite similar. Further understanding of the relaxation mechanism has been obtained from a 
closure look on the heterogeneous dynamics. As figure 3 demonstrates, the rotational non-Gaussian parameter shows 
a dramatic enhancement of hetrogeneous dynamics as the I-N phase boundary is approached. Unlike what is found 
near the gas-liquid critical point [32j, the single-particle dynamics near the I-N phase boundary are observed to be 
strongly affected by the approaching thermodynamic singularity. We have discussed mode coupling theory approaches 
introduced to understand anomalous dynamics observed in this problem. 
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